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Abstract. Wc study equilibrium configurations of swarming biological organisms subject to 
exogenous and pairwise endogenous forces. Beginning with a discrete dynamical model, we derive 
a variational description of the corresponding continuum population density. Equilibrium solutions 
are extrema of an energy functional, and satisfy a Frcdholm integral equation. We find conditions for 
the extrema to be local minimizers, global minimizers, and minimizers with respect to infinitesimal 
Lagrangian displacements of mass. In one spatial dimension, for a variety of exogenous forces, 
endogenous forces, and domain configurations, we find exact analytical expressions for the equilibria. 
These agree closely with numerical simulations of the underlying discrete model. The exact solutions 
provide a sampling of the wide variety of equilibrium configurations possible within our general 
swarm modeling framework. The equilibria typically are compactly supported and may contain 
^-concentrations or jump discontinuities at the edge of the support. We apply our methods to a 
model of locust swarms, which are observed in nature to consist of a concentrated population on the 
ground separated from an airborne group. Our model can reproduce this configuration; quasi-two- 
dimensionality of the model plays a critical role. 

Key words, swarm, equilibrium, aggregation, integrodifferential equation, variational model, 
energy, minimizer, locust 

1. Introduction. Biological aggregations such as fish schools, bird flocks, bacte- 
rial colonies, and insect swarms (21 [2H [30] have characteristic morphologies governed 
by the group members' interactions with each other and with their environment. The 
endogenous interactions, i.e., those between individuals, often involve organisms re- 
acting to each other in an attractive or repulsive manner (TOj l2Z] when they sense 
each other either directly by sound, sight, smell or touch, or indirectly via chemicals, 
vibrations, or other signals. A typical modeling strategy is to treat each individual as 
a moving particle whose velocity is influenced by social (interparticle) attractive and 
repulsive forces [25, 28J. In contrast, the exogenous forces describe an individual's re- 
action to the environment, for instance a response to gravity, wind, a chemical source, 
a light source, a food source, or a predator. The superposition of endogenous and ex- 
ogenous forces can lead to characteristic swarm shapes; these equilibrium solutions 
are the subject of our present study. 

More specifically, our motivation is rooted in our previous modeling study of 
the swarming desert locust Schistocerca gregaria [34.. In some parameter regimes of 
our model (presented momentarily) , locusts self-organize into swarms with a peculiar 
morphology, namely a bubble-like shape containing a dense group of locusts on the 
ground and a flying group of locusts overhead; see Figure [Lltbc). The two are 
separated by an unoccupied gap. With wind, the swarm migrates with a rolling 
motion. Locusts at the front of the swarm fly downwards and land on the ground. 
Locusts on the ground, when overtaken by the flying swarm, take off and rejoin the 
flying group; see Figure [l~T| cd). The presence of an unoccupied gap and the rolling 
motion are found in real locust swarms [TJ |37j . As we will show throughout this 
paper, features of swarms such as dense concentrations and disconnected components 
(that is, the presence of gaps) arise as properties of equilibria in a general model of 
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swarming. 

The model of [34] is 

*i= fc^l^Dpjj -ge z + Ue x , i = l...N, (1.1a) 

q(r) = Ge- r / L -e- r , r y = xj - x*, (1.1b) 

which describes iV interacting locusts with positions Xj . The direction of locust swarm 
migration is strongly correlated with the direction of the wind [3TJ [37] and has little 
macroscopic motion in the transverse direction, so the model is two-dimensional, i.e., 
x s ; = (xi,Zi) where the x coordinate is aligned with the main current of the wind 
and z is a vertical coordinate. As the velocity of each insect is simply a function of 
position, the model neglects inertial forces. This so-called kinematic assumption is 



common in swarming models, and we discuss it further in Section 2.1 



The first term on the right-hand side of (1.1) describes endogenous forces; q{r) 
measures the force that locust j exerts on locust i. The first term of q(r) describes 
attraction, which operates with strength G over a length scale L and is necessary 
for aggregation. The second term is repulsive, and operates more strongly and over 
a shorter length scale in order to prevent collisions. Time and space are scaled so 
that the repulsive strength and length scale are unity. The second term on the right- 



hand side of ( 1.1 1 describes gravity, acting downwards with strength g. The last term 
describes advection of locusts in the direction of the wind with speed U. Furthermore, 
the model assumes a flat impenetrable ground. Since locusts rest and feed while 
grounded, their motion in that state is negligible compared to their motion in the air. 



Thus we add to (1.1 ) the stipulation that grounded locusts whose vertical velocity is 



computed to be negative under (1.1) remain stationary 



As mentioned above, for some parameters, (1.1) forms a bubble-like shape. This 
can occur even in the absence of wind, that is, when [7 = 0; see Figure [OJb). The 
bubble is crucial, for it allows the swarm to roll in the presence of wind. As discussed 
in |34j , states which lack a bubble in the absence of wind do not migrate in the presence 
of wind. Conditions for bubble formation, even in the equilibrium state arising in the 
windless model, have not been determined; we will investigate this problem. 

Some swarming models adopt a discrete approach - as in our locust example 
above - because of the ready connection to biological observations. A further ad- 
vantage is that simulation of discrete systems is straightforward, requiring only the 
integration of ordinary differential equations. However, since biological swarms con- 
tain many individuals, the resulting high-dimensional systems of differential equations 
can be difficult or impossible to analyze. Furthermore, for especially large systems, 
computation, though straightforward, may become a bottleneck. Continuum models 
are more amenable to analysis. One well-studied continuum model is that of [S], a 
partial integrodiffcrcntial equation model for a swarm population density p(x, t) in 
one spatial dimension: 

Pt + (pV) x = 0, V(x)= ( q(x-y)p(y)dy. (1.2) 



The density p obeys a conservation equation, and V is the velocity field, which is 
determined via convolution with the antisymmetric pairwise endogenous force q, the 
one-dimensional analog of a social force like the one in (1.1). 
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Fig. 1.1. Numerical simulation of the two-dimensional locust swarm model of \3$ , given 
by lj !■!]) ■ (a) The simulation begins with a randomly distributed initial state, (b) With gravity but no 
wind (U = 0), the swarm's equilibrium is a bubble-like shape on the ground, consisting of a dense, 
grounded group of locusts and an airborne group. The two are separated by a gap that is void of 
insects, (c) For the full simulation with wind, by the time t = 2, the swarm again coheres into a 
bubble and travels to the right with a rolling motion. Individuals in the back of the swarm take off 
to join the flying group and individuals reaching the front of the flying swarm land on the ground, 
where they remain motionless until taking off again, (d) The trajectory of one individual locust from 
t = to t = 20 demonstrates the periodic landing and takeoff. The parameters in are N = 200, 

G = 0.5, L = 10, g = 1 and U = 1. 



The general model (1.2) displays at least three solution types as identified in [5]. 
Populations may concentrate to a point, reach a finite steady state, or spread. In 23 , 
we identified conditions on the social interaction force q for each behavior to occur. 
These conditions map out a "phase diagram" dividing parameter space into regions 
associated with each behavior. Similar phase diagrams arise in a dynamic particle 
model (T3] and its continuum analog [T2] • Models that break the antisymmetry of q 
(creating an asymmetric response of organisms to each other) display more compli- 
cated phenomena, including traveling swarms [26 . 

Many studies have sought conditions under which the population concentrates to 
a point mass. In a one-dimensional domain, collapse occurs when the force q is finite 
and attractive at short distances [5]. The analogous condition in higher dimensions 
also leads to collapse [3HZ] ■ One may also consider the case when the velocity includes 
an additional term describing an exogenous force, 

V(x)= [ q(x~y)p(y)dy + f(x). (1.3) 
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In this case, equilibrium solutions consisting of sums of point-masses can be linearly 
and nonlinearly stable, even for social forces q that are repulsive at short distances 
[Tol, IT71 ■ These results naturally lead to the question of whether a solution can 
be continued past the time at which a mass concentrates. Early work on a particular 



generalization of (1.2) suggests the answer is yes J2T] [22] . For (1.2 1 itself in arbitrary 
dimension, there is an existence theory beyond the time of concentration |llj . 

Some of the concentration solutions mentioned above are equilibrium solutions. 
However, there may be classical equilibria as well. For most purely attractive q, the 
only classical steady states are constant in space, as shown via a variational formu- 
lation of the steady state problem [5]. However, these solutions are non-biological, 
as they contain infinite mass. There do exist attractive-repulsive q which give rise to 
compactly-supported classical steady states of finite mass. For instance, in simulations 



of (1.2), we found classical steady state solutions consisting of compactly supported 
swarms with jump discontinuities at the edges of the support [23]. In our current 
work, we will find equilibria that contain both classical and nonclassical components. 
Many of the results reviewed above were obtained by exploiting the underlying 



gradient flow structure of (1.3). There exists an energy functional 



W[ P ] = \ [ f p(x)p(y)Q(x-y)dxdy + f F(x)p(x)dx, (1.4) 
£ Jr Jr Jr 

which is minimized under the dynamics. This energy can be interpreted as the con- 
tinuum analog of the summed pairwise energy of the corresponding discrete (particle) 
model [TTJ . We will also exploit this energy to find equilibrium solutions and 
study their stability. 

In this paper, we focus on equilibria of swarms and ask the following questions: 

• What sorts of density distributions do swarming systems make? Are they 
classical or nonclassical? 

• How are the final density distributions reached affected by endogenous inter- 
actions, exogenous forces, boundaries, and the interplay of these? 

• How well can discrete and continuum swarming systems approximate each 
other? 

To answer these questions, we formulate a general mathematical framework for dis- 
crete, interacting swarm members in one spatial dimension, also subject to exogenous 
forces. We then derive an analogous continuum model and use variational methods 
to seek minimizers of its energy. This process involves solution of a Fredholm integral 
equation for the density. For some choices of endogenous forces, we are able to find 
exact solutions. Perhaps surprisingly, they are not always classical. In particular, 
they can involve (5-function concentrations of mass at the domain boundary. 

The rest of this paper is organized as follows. In Section [2] we create the math- 
ematical framework for our study, and derive conditions for a particular density dis- 
tribution to be an equilibrium solution, and to be stable to various classes of pertur- 
bations. 

In Sections [3] andj4] we demonstrate different types of swarm equilibria via ex- 
amples. In Section [3J we focus on purely repulsive endogenous interactions. We 
consider a bounded domain with no exogenous forces, a half-line subject to gravita- 
tional forces, and an unbounded domain subject to a quadratic exogenous potential, 
modeling attraction to a light, chemical, or nutrient source. For all three situations, 
we find exact solutions for swarm equilibria. For the first two examples, these equilib- 
ria consist of a density distribution that is classical in the interior of the domain, but 
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contains (5-functions at the boundaries. For the third example, the equilibrium is com- 
pactly supported with the density dropping discontinuously to zero at the edge of the 
support. For all three examples, we compare analytical solutions from the continuum 
framework to equilibria obtained from numerical simulation of the underlying discrete 
system. The two agree closely even for small numbers of discrete swarm members. 
Section [4] is similar to Section [3j but we now consider the more complicated case of 
endogenous interactions that are repulsive on short length scales and attractive over 
longer ones; such forces are typical for swarming biological organisms. 

In Section [3J we revisit locust swarms, focusing on their bubble-like morphology 
as described above, and on the significance of dimensionality. In a one-dimensional 
model corresponding to a vertical slice of a wide locust swarm under the influence 
of social interactions and gravity, energy minimizers can reproduce concentrations of 
locusts on the ground and a group of locusts above the ground, but there cannot 
be a separation between the two groups. However, a quasi-two-dimensional model 
accounting for the influence of the swarm's horizontal extent does, in contrast, have 
minimizers which qualitatively correspond to the biological bubble-like swarms. 

2. Mathematical formulation. 

2.1. Discrete model. Consider N identical interacting particles (swarm mem- 
bers) in one spatial dimension with positions Xi. Assume that motion is governed 
by Newton's law, so that acceleration is proportional to the sum of the drag and 
motive forces. We will focus on the case where the acceleration is negligible and 
the drag force is proportional to the velocity. This assumption is appropriate when 
drag forces dominate momentum, commonly known in fluid dynamics as the low 
Reynolds number or Stokes flow regime. In the swarming literature, the resulting 
models, which are first-order in time, are known as kinematic. Kinematic models 
have been used in numerous studies of swarming and collective behavior, including 

m ma emb 123 m hi 123 ) . 

We now introduce a general model with both endogenous and exogenous forces, 



as with the locust model (1.1 1. The endogenous forces act between individuals and 
might include social attraction and repulsion; see [28] for a discussion. For simplicity, 
we assume that the endogenous forces act in an additive, pairwise manner. We also 
assume that the forces are symmetric, that is, the force induced by particle i on 
particle j is the opposite of that induced by particle j on particle i. Exogenous forces 
might include gravity, wind, and taxis towards light or nutrients. 
The governing equations take the form 

dx ' 

= Vi(x!,...,x N ), (2.1a) 

at 

N 

Vi(x 1 ,...,x N ) = ^2'mq(xi-x j )+f(xi), m = M/N. (2.1b) 
i=i 

Eventually we will examine the governing equations for a continuum limit of the 
discrete problem. To this end, we have introduced a social mass m which scales the 
strength of the endogenous forces so as to remain bounded for N — > oo. M is the 



total social mass of the ensemble. Eq. (2.1b) defines the velocity rule; mq is the 



endogenous velocity one particle induces on another, and / is the exogenous velocity. 
From our assumption of symmetry of endogenous forces, q is odd and in most realistic 
situations is discontinuous at the origin. 
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Each force, / and q, can be written as the gradient of a potential under the rela- 
tively minor assumption of integrability. As pointed out in [55], most of the specific 
models for mutual interaction forces proposed in the literature satisfy this require- 
ment. Many exogenous forces - including gravity and common forms of chemotaxis 



- do so as well. Under this assumption, we rewrite (2.1 1 as a gradient flow, 
dx 

-^-=V i {x 1 ,...,x N ) = -V 4 W(xi,...,xjv), (2.2) 
where the potential W is 

N N N 

W(asi, . . . ,x N ) = - mQ{x t - Xj ) + F ( x k), (2.3a) 



2 

i=l j=l k=l 



Q(z) = -J q(s)ds, F(z) = -J f(s)ds. (2.3b) 

The double sum describes the endogenous forces and the single sum describes the 
exogenous forces. Also, Q is the mutual interaction potential, which is even, and F is 



the exogenous potential. The flow described by (2.2) will evolve towards minimizers 
of the energy W. 

Up to now, we have defined the problem on K. In order to confine the problem 
to a particular domain Q, one may use the artifice of letting the exogenous potential 
F tend to infinity on the complement of fi. 

While this discrete model is convenient from a modeling and simulation stand- 



point, it is difficult to analyze. Presently, we will derive a continuum analog of (2.1 ). 
This continuum model will allow us to derive equilibrium solutions and determine 
their stability via the calculus of variations and integral equation methods. 

2.2. Continuum model. To derive a continuum model, we begin by describing 
our evolving ensemble of discrete particles with a density function p(x,t) equal to a 
sum of 5-functions. (For brevity, we suppress the t dependence of p in the following 
discussion.) Our approach here is similar to [5]. These S- functions have strength m 
and are located at the positions of the particles: 



i=i 



p(x) =^2m5(x - Xi). (2.4) 

The total mass is 

M = I p(x) dx = mN, (2.5) 



n 



where is the domain of the problem. Using ( |2.4[ ), we write the discrete velocity 
Vi(xi, . . . , xn) in terms of a continuum velocity V(x). That is, we require Vi(x\, . . . , xn) 
V(xi) where 



V(x)= / q(x-y)p(y)dy + f(x). (2.6) 
Jn 

By conservation of mass, the density obeys 



Pt + (pV) x = 0, 



(2.7) 
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with no mass flux at the boundary. We now introduce an energy functional W[p] 
which is analogous to the discrete potential W in (2.3a I: 



W[ P ] = 



p(x)p(y)Q(x - y)dxdy + / F(x)p(x)dx 



Q.JQ. 



This expression follows from the discrete potential, 

W[p] =mW(xx,...,x N ), 



(2.8) 



(2.9) 



remembering the (5-function definition of the density (2.4). The rate of energy dissi- 
pation is 



d -™ = -l /H,>m,>!^, 



(2.10) 



where we assume that there is no mass flux at the boundary of tt. A consequence 
of this boundary condition is that under some conditions, mass may concentrate at 
the boundary of the domain, and we will later see this manifest. Since energy is 
dissipated, we conclude that stable equilibria correspond to minimizers of W. 

is defined for any density distribution p, not 



Imagine now that the energy (2 



just ensembles of (^-functions. We will find minimizers for the continuous energy (2.8) 



and show that they approximate solutions to the discrete problem in the limit of large 
N. To establish a correspondence between the two frameworks, consider a continuous 
distribution p c with total mass M. Define the cumulative density function 



/ Pc(s)ds, 

J Xn 



where the dummy coordinate xq is taken to the left of the support of p c 
discrete approximation of N ^-functions, 



N 



Pd(x) = ^2mS(x - Xi). 



(2.11) 



We seek a 



(2.12) 



i=l 



The associated cumulative density function "F^ is 



x < x\ 

m[l/2 + (£-l)] x = x u i = l,...,N 
n Xi < x < Xi+i, i = 1, . . . , N — 1 



(2.13) 



M 



X > Xjy, 



where we have used the convention that integrating up to a (5-function yields half 
the mass of integrating through it. To establish our correspondence, we require that 
^c{xi) — ^d{xi), which in turn determines the particle positions X{. As N — > oo 
for fixed M, this step function converges uniformly to The correspondence 
goes in the opposite direction as well. We can begin with an ensemble of (^-functions 
pd placed at the positions of discrete swarm members, as shown in Figure 2.1 [&). 



We can find the corresponding cumulative density ^ via (2.13) and interpolate to 
construct the continuum cumulative density ^S c . The functions ^d and ^ c are shown 
as the dotted red step function and the blue curve, respectively, in Figure [2Tjb) . We 
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2 4 6 

Fig. 2.1. Schematic correspondence between discrete and continuum systems in one spa- 
tial dimension, (a) Positions of individual swarm members, modeled as a sum of S -functions 
per l[2.12ty . (b) Cumulative density distributions. The discrete cumulative density 9^(x) appears 
as a step function (dotted red curve) per j[2.1Sty . The continuum cumulative density *D c (x) (solid 
blue curve) is obtained by interpolating ^^(x) such that the two agree at the positions of the swarm 
members (green dots), (c) Continuum density p c (%) corresponding to (a), obtained by differentiating 
the cumulative density ty c (x) in (b). 



may then differentiate to find an approximation p c , as shown in Figure 2.1 ^c). We 



use this correspondence in Sections [3] through [5] to compare analytical results for the 
continuum system (2.7) with numerical simulations of the discrete system (2.1 1. 



We close this subsection by reiterating why we have made a correspondence be- 
tween the discrete and continuum systems. We use the continuous framework to 
find equilibrium solutions analytically via variational and integral equation methods. 
The correspondence above allows a direct comparison to numerical simulation of the 
discrete system. 

2.3. Local and global minimizers. We use a variational calculation to deter- 
mine conditions for a density distribution to be a minimizer of W. Our starting point 



is the energy functional (2.8) subject to the mass constraint (2.5). Let 



p{x) = p + ep(x). 



(2.14) 



Here p is an equilibrium solution of mass M and ep is a small perturbation of zero 
mass, so we have 



p(x) dx = M, 
p(x) dx = 0. 



(2.15a) 
(2.15b) 
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Inspired by biological observations of swarms, we focus on equilibria with finite extent 
and take the support of p to be a finite subset of the domain SI. We refer to the support 
of p as VLp. This assumption, combined with the fact that the density is nonnegative, 
restricts the perturbation p to be nonnegative on Jlp C , the complement of VLp. 
We write 



W[p] = W[p] + eW^p] + e 2 W 2 [p,p}, 



(2.16) 



where W\ and W2 are the first and second variations respectively. This expression is 



exact because W is quadratic in p (see eq. 2.8 1. We analyze these variations to deter- 



mine necessary and sufficient conditions for a candidate solution p to be a minimizer 
of W. 

Our strategy is to consider two classes of perturbations p. First, we consider 
perturbations whose support lies in VLp. In order for p to be extremal, W\ must 



vanish. For it to be a minimizer, W2 must be positive. Since (2.16) is exact, W2 > 



guarantees that p is both a local and global minimum with respect to this first class 
of perturbations. 

The second class of perturbations we consider (of which the first class is a subset) 
consists of perturbations on the entire domain f2. As mentioned above, these pertur- 
bations must be nonnegative in f2p C in order to maintain positivity of p. A necessary 
condition for p to be a local minimizer is that W\ > for this class of perturbations. 
If in addition W2 remains positive for this larger class of perturbations, then p is a 
global minimizer as well. 



We derive the first variation W\ by substituting (2.14) into (2.16) and expanding 
to first order in e, which yields 



Wi[p,p] = I p 
In 



Q(x - y)p{y) dy + F(x) 



dx. 



(2.17) 



Consider the first class of perturbations, whose support lies in Jlp. As the perturbation 
p is arbitrary and of zero mass, for the first variation W\ to vanish, it must be true 
that 



/ Q(x - y)p(y)dy + F(x) = A for 1 e fi p -. 



(2.18) 



The same result can be found by a Lagrange multiplier argument including the con- 
stant mass as a constraint. The multiplier A has a physical interpretation; it is the 
energy per unit mass an additional test mass would feel due to the exogenous potential 
and the interaction with p. 

Thus far, we have shown that a necessary condition for p to be an equilibrium 
solution is that it satisfies the Fredholm integral equation of the first kind for the 
nonnegative density p, 



l[p(x)] = X-F(x), I[p(x)] 



Q(x - y)p{y) dy, 



(2.19) 



as well as the mass constraint (2.15a) 



In order for p to be a minimizer with respect to the first class of perturbations, 

(2.20) 



the second variation must be positive. Substituting (2.14) into (2.16) yields 
W 2 \p,p\= / Q{x-y)p{x)p(y)dxdy. 
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The sign of W% can be assessed in a number of ways. 

We first derive a sufficient condition on the Fourier transform of Q for Wi to be 
nonnegative for the first class of perturbations. Define the Fourier transform 



Q(k)= / Q(x)e 



Then we have 



W 2 [p, p] 



1 

2^ 



dx. 



Q(x - y)p(x)p(y)dxdy, 
Q{x-y)p{x)p{y)dxdy, 
p{x) [Q(x) * p(x)] dx, 
\P(k)\ 2 Q{k) dk. 



Sip •* il p 

OO />OG 



-oo J — oo 

oo 



(2.21) 

(2.22a) 
(2.22b) 
(2.22c) 
(2.22d) 



We have used the fact that p is compactly supported to extend the range of integration 
to infinity. We have also used the convolution theorem, Parseval's theorem, and the 
fact that p is real. We see, then, that Q(k) > is a sufficient condition for W% > 
(assuming a nontrivial perturbation). As shown in |23| . this condition is actually 
equivalent to that for the linear stability of a constant density state in the absence of 
exogenous forces. 

A necessary and sufficient condition for W2 > for the first class of perturbations 
comes from considering the spectrum of I in tip. Note that (2.201 may be written as 



W 2 [p,p] 



l[p]p(x)dx = (p,X[ 



(2.23) 



where the angle brackets denote the usual L 2 inner product on flp. If the eigenvalues 
of the integral operator X are positive, then W2 > (again assuming a nontrivial 
perturbation). 

We now turn to the second class of perturbations, which have support in Q and 
which are positive in Qp C . To analyze these perturbations, we extend the definition 
of the constant A to a function A(x) that is defined over all of fi. We set 



A(x) 



Q{x-y)p{y) dy + F(x). 



(2.24) 



Trivially from (2.181, A(x) = A for x € fip. We now rewrite the first variation as 

Wi[p,p}= / p(x)A(x)dx, 



(2.25) 



directly from (2.17). Remembering that p > in flp C and that p has zero mass in Q, 
we see that a necessary and sufficient condition for W\ > is A(x) > A in £l p c 1 that 
is 

A(x)= I Q{x-y)p{y)dy + F(x)> A fonefl/. (2.26) 



A PRIMER OF SWARM EQUILIBRIA 



11 



Physically, this guarantees that a parcel of mass transported from VLp to its comple- 
ment increases the total energy. 

If we wish for p to be a minimizer with respect to the second class of perturbations, 
it suffices for W2 > (for example, by having Q(k) > 0). However, this condition is 
not necessary. The necessary condition is that W\ [p, p] + W2 [p, p] > for nontrivial 



perturbations, which follows from (2.16) being exact. 

To summarize, we have obtained the following results: 



Equilibrium solutions p satisfy the Fredholm integral equation (|2.19|) and the 
mass constraint ( 2.15a[ ). 



The solution p is a local and global minimizer with respect to the first class 
of perturbations (those with support in Qp) if W 2 in (2.201 is positive. 



The solution p is a local minimizer with respect to the second (more general 



zero-mass) class of perturbations if p satisfies (2.26). If in addition W2 is 



positive for these perturbations, then p is a global minimizer as well. 



In practice, we solve the integral equation (2.19) to find candidate solutions. Then, 
we compute A(x) to determine whether p is a local minimizer. Finally, when possible, 
we show the positivity of W2 to guarantee that p is a global minimizer. 

2.4. Swarm minimizers. As the continuum limit replaces individual particles 
with a density, we need to make sure the continuum problem inherits a physical in- 
terpretation for the underlying problem. If we think about perturbing an equilibrium 
configuration, we note that mass cannot "tunnel" between disjoint components of 
the solution. As such we define the concept of a multi-component swarm equilibrium. 
Suppose the swarm's support can be divided into a set of m disjoint, closed, connected 
components {^i}, that is 



flp = f2i U 2 U • • • U Q m , QinQj = 0, i^j. 



(2.27) 



We define a swarm equilibrium as a configuration in which each individual swarm 
component is in equilibrium, 



Q(x - y)p(y) dy + F(x) = A, for i = l,...,m. (2.28) 
We can still define A(x) in Q 

A(x) = l{p{x)\ + F(x) = [ Q{x- y)p{y) dy + F(x), 



(2.29) 



but now A(x) — Xi in f2j. We can now define a swarm minimizer. We say a swarm 
equilibrium is a swarm minimizer if A(x) > for some neighborhood of each compo- 
nent Vli of the swarm. In practice this means that the swarm is an energy minimizer 
for infinitesimal redistributions of mass in the neighborhood of each component. This 
might also be called a Lagrangian minimizer in the sense that the equilibrium is a 
minimizer with respect to infinitesimal Lagrangian deformations of the distributions. 

It is crucial to note that even if a solution p is a global minimizer, other multi- 
component swarm minimizers may still exist. These solutions are local minimizers and 
consequently a global minimizer may not be a global attractor under the dynamics 

of pTfb. 
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M 
2(1 + d) 



6{x + d) 



M J (x - d) 
2(1 + d) M5{x) 



o 



o 



(a) 



M 
2(1 + d) 



(d) 



P.(x') =0 





Fig. 3.1. Schematic solutions of the minimization problem defined by and \2. 1 5a\ when 

the potential Q(x) is chosen to be the purely repulsive Laplace potential \3. The total mass 
is M. Minimizers p(x) are shown in (a,b,d,e) and two exogenous potentials F(x) are shown in 
(c,f). (a) The minimizer for the bounded domain O = [— d, d] with F(x) = is a constant internal 
density distribution with 8-functions at the boundary. See Section \3.3\ (b) The minimizer for the 
unbounded domain Q = R with an exogenous quadratic potential F(x) of strength 7 is a compactly 
supported, downward-facing parabola section. See Section \ 3.5\ (c) The quadratic potential F{x) 
corresponding to (b). (d) The minimizer for the half-line Q = [0, 00] with an exogenous gravitational 
potential F(x) = gx is a S-function at the origin for sufficiently large gravity g. (e) Like (d), but for 
smaller g, in which case the minimizer is a 5-function at the origin, is linear for a finite interval 
connected to the origin, and drops dis continuously to zero outside of that interval. See Section 
\3.4\ (f) The gravitational potential F(x) corresponding to (d,e). 



3. Examples with a repulsive social force. In this section we discuss the 
minimization problem formulated in Section [2j It is helpful for expository purposes 
to make a concrete choice for the interaction potential Q. As previously mentioned, 
in many physical, chemical, and biological applications, the pairwise potential Q is 
symmetric. Additionally, repulsion dominates at short distances (to prevent collisions) 
and the interaction strength approaches zero for very long distances. A common choice 
for Q is the Morse potential with parameters chosen to describe long-range attraction 
and short-range repulsion [28j . For the remainder of this section, we consider a simpler 
example where Q is the Laplace distribution 

Q{x) = e-W, (3.1) 

which represents repulsion with strength decaying exponentially in space. 
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When there is no exogenous potential, F(x) — 0, and when the domain is infinite, 
e.g., = K, the swarm will spread without bound. The solutions asymptotically 
approach the Barenblatt solution to the porous medium equation as shown in |23j . 
However, when the domain 57 is bounded or when there is a well in the exogenous 
potential, bounded swarms are observed both analytically and numerically, as we 



will show. Figure 3.1 shows solutions p(x) for three cases: a bounded domain with 
no exogenous potential, a gravitational potential on a semi-infinite domain, and a 
quadratic potential well on an infinite domain. In each case, a bounded swarm solution 
is observed but the solutions are not necessarily continuous and can even contain 5- 
function concentrations at the boundaries. 

We discuss these three example cases in detail later in this section. First, we 
will formulate the minimization problem for the case of the Laplace potential. We 
will attempt to solve the problem classically; when the solution has compact support 
contained within the domain we find solutions that are continuous within the support 
and may have jump discontinuities at the boundary of the support. However, when 
the boundary of the support coincides with the boundary of the domain, the classical 
solution may break down and it is necessary to include a distributional component 
in the solution. We also formulate explicit conditions for the solutions to be global 
minimizers. We then apply these results to the three examples mentioned above. 

3.1. Classical solutions to the integral equation. Recall that for p to be 



a steady solution, it must satisfy the integral equation (2.19) subject to the mass 



constraint (2.15a). For p{x) to be a local minimizer, it must also satisfy (2.26) 



A(a:)= / Q(x — y)p(y) dy + F(x) > X. (3.2) 
Finall y, recall that for a solution p to be a global minimizer, the second vari ation 



(2.20) must be positive. We saw that if Q(k) > 0, this is guaranteed. For (3.1), 
Q(k) = 2/(1 + k 2 ) and so for the remainder of this section, we are able to ignore the 
issue of W<l- Any local minimizer that we find will be a global minimizer. 

Additionally, for the remainder of this section, we restrict our attention to cases 
where the support of the solution fip is a single interval in ft; in other words, the 
minimizing solution has a connected support. The reason that we are able to make 
this restriction follows from the notion of swarm minimization, discussed above. In 
fact, we can show that there are no mult i- component swarm minimizers for the Laplace 
potential as long as the exogenous potential F(x) is convex, that is, F"(x) > on 
fl. To see this, assume we have a swarm minimizer with a at least two disjoint 
components. Consider A(x) in the gap between two components so that x £ fip C . We 
differentiate A(x) twice to obtain 

A"(x)= [ Q"{x-y)p{y)dy + F"{x). (3.3) 



Note that Q"(x — y)>0asx^yin VLp C . F"(x) > by assumption. Consequently, 
A"(x) > in f2p C and so A(x) is convex upwards in the gap. Also, A(x) = Xi at the 
endpoints of the gap. We conclude from the convexity that A(x) must be less than Xi 
near one of the endpoints. This violates the condition of swarm minimization from 
the previous section, and hence the solution is not a swarm minimizer. Since swarm 
minimization is a necessary condition for global minimization, we now, as discussed, 
restrict attention to single-component solutions. 
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For concreteness, assume the support of the solution p is Op = [a,/?]. We 



transform the integral equation (2.19) for p into a differential equation by noting 
that Q(x) is the Green's function of the differential operator C = d xx — 1, so that 
CQ(x) — —2S(x), where 6(x) is the Dirac (5-function. Applying C to both sides 
of (2.19) yields the Dirac (5-function under the integral. Integration and the sifting 



property of the 5 function lead to 



p(x)=fk(x) = -~/:[\-F(x)] = ~-~[F-F'Xx)] inn,. (3.4) 



Equation (3.4) is a necessary condition on a solution p(x) but not sufficient. In 



integral equation (2.19). We find 



fact, to verify the candidate solution (3.4), we must substitute back into the governing 



T[p*(x)] = I ( 
In 



X-F(x) + —[F(a)-F'(a)-X] 



(3.5a) 
(3.5b) 



+ -Z-[F(P)+F J (0)-\], 



where the terms other than A — F(x) must vanish for (3.4) to be a solution. These 
error terms are spanned by {e x ,e~ x } which is the null space of C. The constraint 
that these terms must vanish leads to the conditions 



F(a) - F'(a) = A, F(/3) + F'(/3) = A. 



(3.6) 



For a given F(x), both these conditions will only be satisfied for particular choices of 
a and (3. The allowed a and (3 then specify A which in turn determines the mass M 
through (2.15a). As we discuss below, the solution derived here is sometimes, but not 



always, a true minimizer of the energy W . 

3.2. Functional minimization viewpoint and nonclassical solutions. What 



we have presented above is the classical view of the solution of (2.19). The difficulty 



with this calculation is that we expect on physical grounds that for every mass M 
and interval a < x < j3 there should be a minimizer. It is easy to see that if F and Q 
are bounded from below, the energy W is bounded from below also. Let 



Fmin = minF(x), 



min Q(x). 

x>0 



(3.7) 



Then directly from (2.8) 



W > -M 2 Q mm + MF mm . 



(3.8) 



Since W is bounded from below, solutions to the minimization problem exist in the 
space of measure valued functions |39j . Sometimes these solutions are not classical; 
when there is finite attraction at small distances, mass can concentrate in a (5-function 
[SJQT]. For the more biologically relevant case of repulsion at short scales, in free space 
and in the absence of an external potential, solutions are classical for all time [§]. 
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3.2.1. Absence of (^-concentrations interior to a domain. We now show 
that (under sufficient hypotheses) a S- function cannot occur in the interior ofOg. 
Though we have restricted ourselves (above) to considering the Laplace potential (3.1 1, 
this is true more generally. 

In fact, suppose that at a point xq in the interior of f2, 

• the exogenous potential F(x) is twice-differentiable in a neighborhood of x$, 

• Q(x) is twice-differentiable for x > 0, 

• Q{x) has nonzero repulsion at short distances, that is, 

lim q(x) = ±K. (3.9) 

Then a minimizer does not contain a S- function at x — xq. 

To see this, we will compare the energy of a density containing a (5-function at 
x = Xq to one where the (5-function has been replaced by a narrow, unit mass top-hat 
distribution 5 s (x — Xq), 

[0 \x\>e, 

and show that the energy is reduced by 0(e). 

Consider a candidate distribution with a (5-function of mass m at Xq, p(x) — 
Po(x) + mS(x — Xq). Define the change of energy 

AW = W[po + mS s (x - x )] - W[p a + mS(x - x Q )]. (3.11) 



By direct calculation using (2.8), 



AW='-^-l I 5 £ {x-x )Q{x-y)S £ (y-x )dxdy-m 2 Q(0)/2 (3.12) 
Po(y)Q(x - y)[S e (x - x ) - 5(x - x )] dxdy 

IQJU 

+ m I F(x) [S E (x — xo) — 5(x — xq)] dx. 
Jn 

Note that near x = 0, Q(x) = Q(0) - K\x\ + 0(x 2 ) and F(x) = F(x ) + F'(x )(x - 
xo) + 0(x 2 ). Expanding to 0(e 2 ), only the first term of (3.12) persists, yielding 

Km 2 1 f e f E 1 

aw = — ^4^7 J \ x -y\ dxd y + °( £2 ) = -z £m2K + °( £2 )- ( 3 - 13 ) 

For K > 0, corresponding to nonzero repulsion at short distances, AW < 0, indicating 
that the energy is reduced by replacing the (5-function with a narrow top-hat. 

We conclude that (5-functions may not occur in the interior of under the as- 
sumed conditions. However, the above reasoning breaks down for (5-functions on the 
domain boundary, where they may not be replaced by a narrow, symmetric top-hat 
distribution. 

3.2.2. (5-concentrations on a domain boundary. Based on the result above, 
we introduce a candidate solution 



p(x) = p*{x) + A5{x - a) + B6(x - P) in Q, (3.14) 
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where A and B are to be determined. The classical solution is supplemented with 
5-functions at the boundary of Vlp. We will show, in agreement with the calculation 
above, that A and B must vanish unless the boundary of the support, f2p, coincides 
with the boundary of the domain, CI. 



We verify the new candidate solution (3.14) by substituting back into the govern 



ing equation (2.19). We find 

l[p(x)}= [ e-l^l {p{y) + A6(y - a) + B6(y - /?)} dy, 



A - F(x) + ^—[2A + F(a) - F'(a) - A] 



-[2B + F(0)+F'(fl)-\], 



(3.15a) 
(3.15b) 



where the terms other than A — F(x) must vanish for (3.14) to be a solution. This con- 
straint leads to the following conditions on A and B, the coefficients of the nonclassical 
parts of the solution: 



.4 



X — F(a) +F'(a) 



B 



A - F(p) - F'(p) 



(3.16) 



2 ' 2 

Because A is undetermined as of yet, we can find a solution for any mass M . Substi- 



tuting the solution p{x) into the mass constraint (2.15a) yields 



M = p(x) dx, 

= A + B + \ { -P-^-\ f F{x)dx + \[F>{P)-F>{a)], 



AM 



(j3 -a) 



F(a)+F(/3)+ I F(x)dx 



(3.17a) 
(3.17b) 

(3.17c) 



Solving for A in terms of M yields 

2M + F(a) + F(f3) + J Q _ F(x) dx 



A = 



2 + f3-a 



(3.18) 



We've shown that for any Hp = [a, 0] and any mass M we can find a solution 



to (2.19) with p(x) smooth in the interior and with a concentration at the endpoints. 



However, we haven't yet addressed the issue of p(x) being non-negative, nor have we 
considered whether it is a minimizer. 

We next consider whether the extremal solution p is a minimizer, which involves 
the study of ( 3.2 ) . We present a differential operator method that allows us to compute 



A(x) and deduce sufficient conditions for p to be a minimizer. 

We start by factoring the differential operator C = d xx — 1 = T) + 'D^ = r D~T> + 
where T) ± = d x ± 1. Applying these operators to the interaction potential Q, we see 
that 



V + Q{x-y) = 



2e x ~ v y > x 







y <x, 



V-Q(x-y) = 







y>x 
-2e v ~ x y < x. 



(3.19) 
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Substituting p in (3.14) into our definition of A(x) in (3.2) yields 



Q(x - y)p(y) dy + AQ(x - a) + BQ(x -(3)+ F(x) = A{x). 



Now consider applying V to (3.20) at a point x in Sip. We see that 



p(y) dy - 2Ae a ~* + V~[F(x)] = T>~ [A] , 



(3.20) 



(3.21) 



where we've used the fact that A(x) = A in Sip. If we let x 
to zero, the integral term vanishes and 

- 2A + F'(a) - F(a) = -A. 



and let z decrease 



(3.22) 



Solving for A yields the first half of (3.16). A similar argument near x — f3 yields the 
value of B. 

Assuming a does not coincide with an endpoint of SI, we now consider the region 
x < a, which is to the left of the support Sip. Again, applying T>~ to (2.19) simplifies 



the equation; we can check that both the integral term and the contribution from the 
^-functions are annihilated by this operator, from which we deduce that 



V~[F(x) - A(x)} = 



F(x) - A(x) = Ce x 



(3.23) 



where C is an unknown constant. A quick check shows that if F(x) is continuous, then 
A(x) is continuous at the endpoints of Sip so that A(a) = A. This in turn determines 
C, yielding 



A{x) = F(x) + [A - F(a)}e x - a 
A similar argument near x — (3 yields 

A(x) = F{x) + [A - F(/3)]e /3 -* 



for x < a. 



for x > p. 



(3.24) 



(3.25) 



As discussed in Section [2. 3| for p{x) to be a minimizer we wish for A(a;) > A for x > (3 
and x < a. A little algebra shows that this is equivalent to 



F(x)e~ x - F{a)e- 1 

e -x _ e -a 

F(x)e x - F{P)et 



> A for x < a, 



> A for x > (3. 



(3.26a) 
(3.26b) 



If a and f3 are both strictly inside SI, then (3.26) constitutes sufficient conditions for 



the extremal solution p to be a global minimizer (recalling that W2 > 0). We may 



also derive a necessary condition at the endpoints of the support from (3.26). As x 



increases to a, we may apply L'Hopital's rule and this equation becomes equivalent to 
the condition A < 0, as expected. A similar calculation letting x decrease to /? implies 
that B < 0. However, since p is a density, we are looking for positive solutions. Hence, 
either A = or a coincides with the left endpoint of SI. Similarly, either B = or 



(3 coincides with the right edge of SI. This is consistent with the result (3.13) which 
showed that S- functions cannot occur in the interior of SI. 
In summary, we come to two conclusions: 

• A globally minimizing solution p contains a (5-function only if a boundary of 
the support of the solution coincides with a boundary of the domain. 

• A globally minimizing solution p must satisfy ( 3.26[ ). 
We now consider three concrete examples for SI and F(x). 
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3.3. Example: Bounded interval. We model a one-dimensional biological 
swarm with repulsive social interactions described by the Laplace potential. We begin 
with the simplest possible case, namely no exogenous potential, F(x) = and a finite 
domain which for convenience we take to be the symmetric interval ft = [—d,d\. 
As F"(x) = we know from (3.3| that the minimizing solution has a connected 
support, i.e., it is a single component. We will see that the minimizing solution has 
an equipartition of mass between ^-functions at the boundaries of the domain and a 



constant solution in the interior, as shown schematically in Figure 3.1 a) 



We now proceed with calculating the solution. From (3.4), we find that 



p*{x) = -. 



The nonclassical solution is 



(3.27) 



p{x) = + AS{x + d) + BS(x - d) in ftp. 



Eq. (3.18) gives A as 



Eq. (3.16) specifies A and B as 



A = B 



M 
l + d' 



M 



2(1 + d) 



(3.28) 



(3.29) 



(3.30) 



Since SI — ftp, it follows that A(x) = A and Wi vanishes according to (2.25). There- 
fore, the solution is a global minimizer. 

The solution p(x) is shown schematically in Figure |3T[ a) 



Figure |3T2[ a) compares 
analytical and numerical results for an example case where we take the total mass to 
be M = 1 and the finite domain to be SI — [—d,d] with d = 1. Cross-hatched boxes 
indicate the boundary of the domain. The solid line is the classical solution p. Dots 



correspond to the numerically-obtained equilibrium of the discrete system (2.1) with 
N = 40 swarm members. The density at each Lagrangian grid point is estimated 
using the correspondence discussed in Section (2.2) and pictured in Figure 2.1 Each 



"lollipop" at the domain boundary corresponds to a (5-function of mass 10/N ■ M = 
1/4 in the analytical solution, and simultaneously to a superposition of 10 swarm 
members in the numerical simulation. Hence, we see excellent agreement between 
the continuum minimizer and the numerical equilibrium even for this relatively small 
number N = 40 of Lagrangian points. 

3.4. Example: Gravitational potential on the half-line. We now consider 
repulsive social interactions and an exogenous gravitational potential. The spatial 
coordinate x > describes the elevation above ground. Consequently, SI is the semi- 
infinite interval < x < oo. Then F(x) = gx with g > 0, shown in Figure [3~T| (f). As 
F"(x) = we know from (3.3) that the minimizing solution has a connected support, 
i.e., it is a single component. Moreover, translating this component downward de- 
creases the exogenous energy while leaving the endogenous energy unchanged. Thus, 
the support of the solution must be ftp = [0,/3], potentially with j3 = 0. In fact, we 
will see that there are two possible solution types depending on g. For strong enough 
values of g, f3 indeed equals zero, and the mass accumulates on the ground x = 0, as 
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1.5 



1.4 



1.3 



(c) 




-0.4 



-0.2 



0.2 



0.4 



Fig. 3.2. Comparison of numerical and exact solutions of the minimization problem \2. 1 9[ ) and 
i 2.1 5a[ ) when the potential Q(x) is the repulsive Laplace potential \3.1\ l. The total mass is M = 1. 
Dots correspond to equilibria of numerical simulations of the underlying discrete system \2. j| ) with 
N = 40 swarm members. The density at the location of each swarm member is estimated using 
the correspondence of Section \2.2\ visualized in Figure \2.1\ The solid curves represent exact solu- 
tions, (a) Bounded domain £7 = [—1,1] with no exogeneous potential, i.e., F(x) = 0. Cross-hatched 
boxes indicate the domain boundaries. The "lollipops" at each end are drawn at an arbitrary height 
and represent 5 -functions each of mass 10/40 in the exact continuum solution, and a superposition 
of 10 swarm members each in the simulation. See Section \ 3. 3\ and Figure \3.1\ a). (b) The half-line 
f2 = [0, oo] with an exogenous gravitational potential F{x) = gx, g = 0.5. The S-function at the 
origin has mass 28/40. See Section \3.4\ and Figure \3.1[ e,f). (c) Unbounded domain Q = R with an 
exogenous quadratic potential F(x) = -yx 2 , 7 = 1. The compactly supported solution drops discon- 
tinuously t o zer o, and the vertical axis is broken for convenience of visual display. See Section \3.5\ 
and Figure\3~l]fbc). 
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shown schematically in Figure |3.l[ d). For weaker g, the mass is partitioned between a 
J-function on the ground and a classical solution for x > 0, as shown in Figure [3~T| e). 
We now proceed with calculating the solution. From (3.4), we find that 

A - gx 



The nonclassical solution is 

p(x) = p* 



p*{x) 



AS(x) + B6(x - 



in Qp. 



Eq. (3.18) gives A as 



Eq. (3.16) specifies A and B as 
A M 



2 + ft 



A = 



2M 

2 + f3 



B 



M 
2 + (3 



|(2 + /3). 



(3.31) 



(3.32) 



(3.33) 



(3.34) 



For p to be a global minimizer, it must be true that B = as shown in Section |3.2 
Solving, we find that (3 = 2-y/M/g - 2. 

For M > g, we see that ft > in which case the minimizing solution is 



p(x) 







IgM- § (1 + x) + s/gM8{x) 0<x< 2^jMjg-2 

x > 2^/M/g - 2. 



(3.35) 



It follows directly that p(x) > 0. Figure |3.l| (e) shows a schematic drawing of this 
solution. 

We still must consider the condition (3.26b). To see it is satisfied, first note that 
A = g(ft + 1). Then ( |3~26b| becomes 

F(x)e x - F(ft)e p 



>9(1 



(3.36) 



We note that by l'Hopital's Rule that equality is obtained as x — > ft. To show the 
inequality holds for x > ft let x = In s and ft = In r (so s > r > 0) . The inequality 
becomes 



sF(lns) — rF(hir) 



> g(l + ft) for s > r. 



(3.37) 



We can interpret the left-hand side as the slope of a chord connecting the points 
(r,rF(lnr)) and (s, sF(his)). Consequently, if the function sF(lns) is concave up- 
ward, the slope of the chord will be increasing as s increases away from r, and the 
inequality will hold. Recalling that F(x) — gx, we compute [sF(lns)]" = g/s which 
is positive, and hence the solution is globally stable. 

For M < g our previous calculation naively implies ft < 0. Since ft cannot be 
negative, the minimizer in this case is a 6- function at the origin, namely p(x) = M8{x), 
shown in Figure 3.1 d) . In this case, A = M from (3.33 ) and A(cc) = Me~ x + gx from 
( |3.25[ ). It follows that 

A(x) = Me~ x + gx >M(l-x) + gx 



M + (g- M)x > M for x > 0. (3.38) 
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The first inequality follows from a Taylor expansion. The second follows from our 
assumption M < g. Since A(x) > A the solution is a global minimizer. 

In summary, there are two cases. When M < g, the globally stable minimizer is a 
(5-function at the origin. When M > g there is a globally stable minimizer consisting 
of a 5-function at the origin which is the left-hand cndpoint of a compactly-supported 
classical swarm. The two cases are shown schematically in Figures |3.1[ de) . Figure 
|3.2[ b) compares analytical and numerical results for the latter (M > g) case with 
M = 1 and g — 0.5. We use N = 40 swarm members for the numerical simulation. 
The numerical (dots) and analytical (line) p* agree, as does the nonclassical part 
of the solution, pictured as the "lollipop" which represents a superposition of 28 
swarm members in the numerical simulation having total mass 28/ N ■ M = 0.7, and 
simultaneously a (5-function of mass \JgM = l/v2 « 0.707 in the analytical solution. 

3.5. Example: Quadratic potential on the infinite line. We now consider 
the infinite domain Q = M with a quadratic exogenous potential well, pictured in 



Figure 3.1 c). This choice of a quadratic well is representative of a generic potential 



minimum, as might occur due to a chemoattractant, food source, or light source. Thus 
F(x) = 7a: 2 where 7 > controls the strength of the potential. As F"(x) = 2~f > 



we know from (3.3) that the minimizing solution is a single component. We take the 



support of the solution to be Clp = [a, 0\. We will find that the compactly supported 
density is classical and has an inverted parabolic shape, shown in Figure [3~T| b). 

Our calculation proceeds as follows. We know from Section |3.2| that because a 



and (3 are assumed to be finite, A = B — and so the solution is classical. From (3.4 1 
we find that 



p*(x) = 



A + 27 - 7a; 2 



(3.39) 



and, from (3.16 1 



A = A - 7« 2 + 2 7 Q = Q B = A - 7 /3 2 - 2 7 /3 = Q 



(3.40) 



Eliminating A from these equations and recalling that /3 > a, It follows that a = —0, 
Hence, the solution is symmetric around the center of the potential. For convenience, 
we now define ft = —a = H. 



Eq. (3.18) gives A in terms of the mass M as 

M + 7 (H 2 + H 3 /3) 



A = 



1 + H 



(3.41) 



However, from the second half of (3.40) we know that 

\ = l{H + 1) 2 -7. 



(3.42) 



Equating these two expressions for A yields that 



V2 7 



1/3 



1. 



(3.43) 



Note that H increases monotonically from with increasing M. 
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Finally, writing the solution p in terms of M, we have 



p(x) 




1-x 2 



x < 



27 



1/3 
1 -1 



N>(W + 1 ) 



1/3 



(3.44) 



1. 



It follows directly that p(x) > 0, and that there is a discontinuity at the edge of the 
support. 

We now must show that the candidate p is a global minimizer, which is done 



3.4 



by dem onst rating (3.26). By the same argumentation used for the right endpoint in 
Section 



it suffices to show that sF(lns) is concave up for s > e H which follows 
directly from algebra. Hence, the solution is a global minimizer. 

The solution p is shown schematically in Figure |3.1| b) . Figure |3.2[ b) compares 
analytical and numerical results for the case with M = 1 and 7=1, with N = 40 
swarm members used for the numerical solution. The numerical (dots) and analytical 
(line) p agree closely. 

4. Examples with a Morse-type social force. In many physical, chemical, 
and biological applications, the pairwise potential Q is isotropic, with a repulsion 
dominating at short distance and the interaction strength approaching zero for very 
long distances. A common choice for Q is the Morse potential 



Q(x) = -GLe-W' L + e~W, 



(4.1) 



with appropriately chosen parameters. Here, x is the distance between particles, 
L > 1 is the characteristic length scale of attraction, and G < 1 is the characteristic 
velocity induced by attraction. We have scaled the characteristic repulsive strength 
and length scale to be unity. In this section we are concerned with the solution of 



(2.15a) and (2.191 with Q given by (4.1) 



The Morse potential has been studied extensively and has become a canonical 
model for attractive-repulsive interactions [14, 24, 28, [34]. A key characteristic of the 



potential (4.1 ) is whether or not the parameters G, L are chosen to be in the H-stable 
or the catastrophic regime; see [33] for a review. Consider (2.1 ) with F(x) = 0. If the 
parameters G, L are chosen in the H-stable regime GL 2 < 1 , then as the number of 
particles N increases, the density distribution of particles approaches a constant, as 
does the energy per particle. Stated differently, the particles form a crystalline lattice 
where the nearest-neighbor distance is approximately equal for all particles (excluding 
edge effects). As more individuals are added to the group, the inter-organism spacing 
is preserved and the group grows to cover a larger spatial region. If the parameters 
are chosen outside of the H-stable regime, i.e., GL 2 > 1, the system is catastrophic. 
In this case, the energy per particle is unbounded as N — » 00, and particles pack 
together more and more closely as N increases. 

Our work in [23j classifies the asymptotic behaviors of the continuum prob- 



lem (2.7) with Morse-type interactions in the absence of an external potential, F(x) = 
0, and on the real line. In the H-stable regime, the continuum model displays spread- 
ing density profiles, while in the catastrophic regime, it forms compactly-supported 
steady states. 

The distinction between catastrophic and H-stable is related to the Fourier trans- 
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form of the potential Q(x). Note that (4.1 ) has Fourier transform 

2 2GL 2 



Q{k) = 



1 + k 2 1 - 
2[(1-GL 2 ) 



(kL) 2 
f L 2 (l 



G)k 2 ] 



(l + k 2 )(l + (kL) 2 ) 



(4.2a) 
(4.2b) 



The condition for H-stability is equiva lent to Q(0) 
Q{k) > for all fc, and from Section 



2.3 



2(1 - GL 2 ) > 0. In this case, 
this is sufficient for Wi > 0, or equiva- 
lently, linear stability of constant density states. In the catastrophic case, Q(0) < 0. 
Intuitively, in the catastrophic case, the constant density state is unstable to long 
wave perturbations. The system is attracted to states of finite extent shorter than 
the length scale of the instability. In the H-stable case, the constant density state 
is stable to perturbations, and initial density profiles spread evenly to become flat. 
We will now study minimizers for the case of Morse-type interactions, and will see 
qualitatively different solutions for the catastrophic and H-stable cases. 

4.1. Classical and nonclassical solutions. We follow the procedure used in 
Section[3j namely we first look for a classical solution on the interior of Q,p = [a, (3] and 
then allow for ^-functions at the boundaries. Once again, we will see that minimizers 
contain <5-functions only when the boundary of the support, flp, coincides with the 
boundary of the domain, Q. 

For convenience, define the differential operators L\ = d xx — l and L 2 = L 2 d xx — 1, 



and apply —C1C2 to (2.19) to obtain 



where 



Thus, we guess 



v{p*)xx + ep* = -CiC 2 [\ - F(x)], xeQp, 



2L 2 (1-G)>0, e = 2(GL 2 -l). 



p(x) = p*(x) + A5(x - a) + B8(x - /3). 



(4.3) 



(4.4) 



(4.5) 



The full solution to the problem is obtained by substituting (4.5) into (2.15a) and 



(2_19| which yields 

P 



Q(x - y)p*(y) dy + AQ(x - a) + BQ{x - 0) = A(x) - F(x), (4.6) 



where A(x) — A in Hp. We begin by considering the amplitudes A and B of the 
distributional component of the solution. We factor the differential operators C± = 
d xx -l= V+V- = V-V+ where V ± = d x ± 1 and C 2 = L 2 d xx -l = Q+ Q~ = Q~ Q+ 
where Q ± = Ld x ± 1. Note that 

V-Q-Q(x-y) = 2L(G-l)6(x-y)-2(L+l) \Ge {y - x)/L - e y ~ x ] H(x-y), (4.7) 



where H is the Heaviside function. Now we apply V Q to (4.6) at a point x in fi^, 
which yields 

2L(G - (a:) - 2{L + 1) f \{Ge^' L - e""*] p* (y) dy (4.8) 

J a 

+ AV-Q-Q(x-a)=V-Q-{\-F(x)}. 
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Taking the limit x — > a + yields 



2L(G - l)p,(a) + 2A(L + 1)(1 - G) = A - V-QrF{x) 



(4.9) 



where we have used the fact that V Q A = A. A similar calculation using the 
operators V + Q + and focusing near x = (3 yields that 



2L(G - l)p*03) + 2B(L + 1)(1 - G) = A - P+g+F(x) 



x=/3 



(4.10) 



Eqs. (4.9) and (4.10) relate the amplitudes of the S- functions at the boundaries to 



the value of the classical solution there. Further solution of the problem requires 
F{x) to be specified. 

In the case where Q, = f2g, solving (4.3) for p*(x) and solving (4.9) and (4.10) 



for A and B yields an equilibrium solution. One must check that the solution is 
non-negative and then consider the solutions stability to determine if it is a local or 
global minimizer. In the case where Qp is contained in the interior of 0, we know 
that A = B = as discussed in Section f3. 2. II We consider this case below. 

Suppose flp is contained in the interior of Q. Then A = B = 0. Following 
Section 3.2 we try to determine when A(x) > A in 17 and when W2 > 0, which 
constitute necessary and sufficient conditions for p to be a global minimizer. We 



apply V~ Or to (4.6) at a point x < a. The integral term and the terms arising from 
the 8 functions vanish. The equation is simply V~ Q~ {A(x) — F(x)} = 0. We write 
the solution as 



A(x) = F(x) + k ie x - a + k 2 e {x ~ a)/L . 



(4.11) 



The two constants fci^ are determined as follows. From (4.6), A(x) is a continuous 
function, and thus 



A(a) = F(a) + k\ + k 2 = A. 



(4.12) 



We derive a jump condition on the derivative to get another equation for k\^- We 
differentiate (|4.6[) and determine that A'(x) is continuous. However, since A(x) = A 



for x € fip, A' (a) = 0. Substituting this result into the derivative of (4.11 ) and letting 
x increase to a~ wc find 



A' (a") = F'(a) + h + k 2 /L = 0. 



(4.13) 



The solution to (J4T2J) and ( |4.13[ ) is 

h 1 



L 



-{F{a)-LF'{a)-X}, 



fc 2 = 



L 



L - 



-{-F{a)+F'{a) + X}. 



(4.14a) 
(4.14b) 



Now that A(x) is known near x — a we can compute when A(x) > A, at least 
near the left side of £lp. Taylor expanding A(x) around x — a, we find 



A(x) « A + ~{P-Q-F{x)\ x=a - A}(x - af + 



(4.15) 
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The quadratic term in (4.15) has coefficient 



V-Q-F(x)\ 
2L 

(1-G)p.(a) >0, 



(4.16a) 
(4.16b) 



where the second line comes from substituting (4.9) with A = and noting that 



the classical part of the solution p*(a) must be nonnegative since it is a density. 
Furthermore, since we expect p*(a) > (this can be shown a posteriori), we have 
that the quadratic term in ( 4.15| ) is positive. 

A similar analysis holds near the boundary x = 0. Therefore, A(x) > A for x 



in a neighborhood outside of ftp. Stated differently, the solution (4.5) is a swarm 
minimizer, that is, it is stable with respect to infinitesimal redistributions of mass. 



The domain flp is determined through the relations (4.9) and (4.10), which, when 
A = B = 0, become 



2L(G - l)p*(a) = A - V~Q~F{x) 
2L(G - l)p*0S) = A - V+Q+F(x) 



(4.17a) 
(4.17b) 



In the following subsections, we will consider the solution of the continuum system 

0. We consider two cases for 
I, for which 



(2.19) and (2.15a) with no external potential, F{x) 
the Morse interaction potential (4.1): first, the catastrophic case on ft 



the above calculation applies, and second, for the H-stablc case on a finite domain, in 
which case A, B ^ and there are ^-concentrations at the boundary. Exact solutions 
for cases with an exogenous potential, F(x) ^ can be straightforwardly derived, 
though the algebra is even more cumbersome and the results unenlightening. 

4.2. Example: Catastrophic interactions in free space. In this case, F(x) — 
in &2A9\ and Gl? > 1 in (O) so that e > in The solution to (Ob is 



p(x) — p*(x) = C cos(p,x) + Dsiii(px) — A/e, 



where 




GL 2 - 1 
L 2 {l-G)' 



(4.18) 



(4.19) 



In the absence of an external potential, the solution is translationally invariant. Con- 
sequently, we may choose the support to be an interval flp — [-H, H] which is sym- 



metric around the origin. Hence, by symmetry, D — 0. While the solution (4.18) 
satisfies the ordinary differential equation (4.3), substituting into the integral equa- 



tion (2.19) and the mass constraint (2.15a) will determine the constants C, A and 



H. The integral operator produces modes spanned by {coshz, coshz/i}. From these 
modes follow two homogeneous equations for C and A which simplify to 



l + /i 2 



1 + p 2 L 2 



[cos(pH) — psin(p,H)] C 
[cos(pH) — fiL sin(pH)] C 



GL 2 - 1 



GL 2 - 1 



A = 0, 



A = 0. 



(4.20a) 
(4.20b) 
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(a) 



0.15 



0.265 




0.26 



Fig. 4.1. Comparison of numerical and exact solutions of the minimization problem 
and \2. 1 5a| ) when the potential Q(x) is the Morse potential The total mass is M = 1. Dots 

correspond to equilibria of numerical simulations of the underlying discrete system \Z. 1\ with N = 40 
swarm members. The solid curves represent exact solutions, (a) Unbounded domain Q = R with 
parameters F = 0.5, L = 2 in \4- 1\ , which is in the catastrophic parameter regime. See Section 
\4-S\ (b) Bounded domain f2 = [— 1,1] with F = 0.1, L = 2, which is in the H-stable parameter 
regime. The 8-functions at the boundary each have mass 10/40 and the vertical axis is broken for 
convenience of visual display. See Section\4-3\ 



For these equations to have a nontrivial solution for C and A the determinant of the 
coefficient matrix must vanish, which yields a condition specifying H , 



cot(nH) = 



GL - 1 



V(1-G)(GL 2 -1)' 



(4.21) 



The mass constraint (2.15a) yields 



(4.22) 
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Solving (4.20a) and (4.22) for C and A yields the full solution for the coefficients in 



(4.18) and the half-width H of the solution, 



H 



C = 



cot 



[O.tt] 

71/ 



GL — 1 
v/(l-G)(GL 2 

Vg(l 2 -i) 



X 



2(H- 
M(l 



•L + l) 
- GL 2 ) 



L{1 - G) 



(4.23a) 

(4.23b) 
(4.23c) 



H + L + l 

As we have shown in Section [4. 1[ this solution is a swarm minimizer. In fact, the 



solution is also a local minimizer. To see this, note that from (4.11) and (4.14) that 



A(x) 



Ixl < H 



(4.24) 



Since A < 0, we see that A(x) > A for |x| > H, ensuring that the solution is a local 
minimizer. 

While we sus pect that is a global minimizer, this is not immediately apparent 
because Q(k) in (4.2) has mixed sign in this catastrophic case, and hence W2 is of 



indeterminate sign. To establish that is a global minimizer one might study the 
quantity W\ + W2 but we leave this analysis as an open problem. 

Figure [4~T| a) compares analytical and numerical results for an example case with 
total mass M = 1 and interaction potential parameters F = 0.5 and L = 2. The 
solid line is the compactly supported analytical solution p. Dots correspond to the 



numerically-obtained equilibrium of the discrete system (2.1) with N — 40 swarm 
members. 

4.3. Example: H-stable interactions on a bounded domain. As described 
in [23], the asymptotic behavior of (2.7) for the H-stable case is a spreading self- 



similar solution that approaches the well-known Barenblatt solution of the porous 
medium equation. Hence, there is no equilibrium solution for the H-stable case on 
an unbounded domain (one can verify this by considering the analogous problem to 
that of the previous section and showing explicitly that there is no solution). Here, 
we assume a bo unded domain f2 = \— d, d}. As before, F(x) — in (2.19) but now 
GL 2 < 1 in (|4.l|) so that e < in (|4.3|). The classical solution to (|4.3|) is 



p* (x) = C cosh(//x) + D sinh(/ix) — A/e, 



where 



V = 




(4.25) 



(4.26) 



We will again invoke symmetry to assume D — 0. The minimizer will be the classical 
solution together with <5-functions on the boundary, 

p(x) = p*(x) + A5{x + d) + B6(x - d). (4.27) 

Again by symmetry, B — A. Consequently, the solution can be written as 

p = Ccosh(/2x) - A/e + A[6(x + d) + 5(x - d)}. (4.28) 
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Substituting into the integral equation (2.19) and the mass constraint (2.15a! will 
determine the constants C, A and A. The integral operator produces modes spanned 
by {cosh z, cosh z/L}. This produces two homogeneous, linear equations for C, A and 
A. The mass constraint (2.15a I produces an inhomogeneous one, namely an equation 
linear in C, A, and A for the mass. We have the three dimensional linear system 

p,d p-fid 



2 yi-p, i + ji 



L 



- /7 d 



1 — flL 1 + jlL 
2 sinh(/id) 
A 



2(1 - GL 2 ) 
L 

2(1 - GL 2 ) 
d 

1 - GL 2 



'A' 




"0" 


C 







A 




M 



(4.29) 



The solution is 

A 
C 
A 



fi 2 LM [(1 + £)(! + jlL)efi d - (1 - /*)(! - /L^e"^] 



2$ 



AM(1 - A )(1 - A £ ) 
Mfi{\ - GL 2 ) [(1 + AK1 + P-L) + (1 - A)(i - A^)] 



(4.30a) 
(4.30b) 
(4.30c) 



where for convenience we have defined 



$ = (l + A)(l + /ii)(A£ + Ad+A-lK d -(l-il){l-ilL){pL + fld+il + l)e-^ a . (4.31) 

For this H-stable case, Q(k) > which ensures that W2 > for nontrivial per- 
turbations. This guarantees that the solution above is a global minimizcr. 

In the limit of large domain size d, the analytical solution simplifies substantially. 
To leading order, the expressions (4.30) become 



A 



C 



X = 



JxLM 
2d ' 

e-f* d M(l - A)(l - flL) 
d : 

M{1 - GL 2 ) 



(4.32a) 
(4.32b) 
(4.32c) 



Note that C cosh(Ax) is exponentially small except in a boundary layer near each edge 
of f2, and therefore the solution is nearly constant in the interior of Q. 

Figure |4~T| b) compares analytical and numerical results for an example case with 
a relatively small value of d. We take total mass M = 1 and set the domain half- 
width to be d = 1. The interaction potential parameters F = 0.5 and L = 2. The 
solid line is the classical solution p*. Dots correspond to the numerically-obtained 
equilibrium of the discrete system (2.1 ) with N = 40 swarm members. Each "lollipop" 
at the domain boundary corresponds to a (5-function of mass 10/ AT • M — 1/4 in the 
analytical solution, and simultaneously to a superposition of 10 swarm members in 
the numerical simulation. 
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5. Modeling a locust swarm: Examples with a gravitational potential. 

We now return to the locust swarm model of [33] , discussed also in Section [T] Recall 
that locust swarms are observed to have a concentration of individuals on the ground, a 
gap or "bubble" where the density of individuals is near zero, and a sharply delineated 



swarm of flying individuals. This behavior is reproduced in the model (1.1); see 
Figure |D|b). In fact, Figure |l.l[ c) shows that the bubble is present even when the 
wind in the model is turned off, and only endogenous interactions and gravity are 
present. 

To better understand the structure of the swarm, we consider the analogous con- 
tinuum problem. To further simplify the model, we note that the vertical structure 
of the swarm appears to depend only weakly on the horizontal direction, and thus 
we will construct a quasi-two-dimensional model in which the horizontal structure is 
assumed uniform. 

In particular, we will make a comparison between a one-dimensional and a quasi- 
two-dimensional model. Both models take the form of the energy minimization prob- 



lem (2.19 1 on a semi-infinite domain, with an exogenous potential F(x) = gx describ- 
ing gravity. The models differ in the choice of the endogenous potential Q, which 
is chosen to describe either one-dimensional or quasi-two-dimensional repulsion. The 
one-dimensional model is precisely that which we considered in Section [3. 4| There we 
saw that minimizers of the one-dimensional model can reproduce the concentrations of 
locusts on the ground and a group of individuals above the ground, but there cannot 
be a separation between the grounded and airborne groups. We will show below that 
for the quasi-two-dimensional model, this is not the case, and indeed, some minimizers 
have a gap between the two groups. 

As mentioned, the one-dimensional and quasi-two-dimensional models incorporate 
only endogenous repulsion. However, the behavior we describe herein does not change 
for the more biologically realistic situation when attraction is present. We consider 
the repulsion-only case in order to seek the minimal mechanism responsible for the 
appearance of the gap. 

5.1. The quasi-two-dimensional Laplace potential. We consider a swarm 
in two dimensions, with spatial coordinate x = (xi,x 2 )- We will eventually confine 
the vertical coordinate x\ to be nonnegative, since it describes the elevation above 
the ground at x\ = 0. We assume the swarm to be uniform in the horizontal direction 
x 2 , so that p(xi,x 2 ) = p{x\). 

We construct a quasi-two-dimensional interaction potential, 



Q2d{\xx-Vi\) = - Q(\x-y\)dy 2 . (5.1) 

J — OO 

Letting z — x\ — yi and s — xi — 7/2? this yields 

pOO 

Q2d(z) = / Q(v / s 2 +z 2 ) ds. (5.2) 



It is straightforward to show that the two-dimensional energy per unit horizontal 
length is given by 

W 2 d[p\ = \ [ I p{xi)p{yi)Q2D(xi-yi)dx 1 dy 1 + I F(x 1 )p(xi) dx u (5.3) 
* Jq Jq. Jn 

where the exogenous force is F(xi) = gx\ and the domain f2 is the half-line x\ > 0. 



This is exactly analogous to the one-dimensional problem (2.8), but with particles 
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interacting according to the quasi-two-dimensional endogenous potential. Similarly, 
the corresponding dynamical equations are simply ( 2.6 ) and (2.7| but with endogenous 
force q 2D = ~Q' 2 d- 

For the Laplace potential (3.1), the quasi-two-dimensional potential is 



ho(z) = 



ds. 



This integral can be manipulated for ease of calculation, 



Q2d(z) 



-^+^ ds, 



2 

2\z 

2\z\e-^ I + 



e -MVTW dWj 

e -\A u 



yV + 2u 

,tt/2 

2\z\ / e'^ scce sec6 d0, 
Jo 



du. 



(5.4) 

(5.5a) 
(5.5b) 
(5.5c) 
(5.5d) 
(5.5e) 



where (5.5b) comes from symmetry, (|5.5c ) comes from letting s = \z\w, (5.5d) comes 
from letting y/1 + w 2 = u+ 1, and (5.5e) comes from the trigonometric substitution 
to = tan 9. 

From an asymptotic expansion of (5.5d), we find that for small \z\, 
Q 2D (z) pa 2 + (7 - 1/2 - In 2) z 1 + z 2 \n\z\ + O (z 4 In |z|) , 
whereas for large \z\, 



(5.6) 



Q 2D (z) Pd y^\z] e-^ 



1 



15 



G(l/|z| 3 ) 



(5.7) 



\z\ 128z 2 

In our numerical study, it is important to have an efficient method of computing values 
of Qw- In practice, we use (5.6) for small |z|, (5.7) for large \z\, and for intermediate 
values of \z\ we interpolate from a lookup table pre-computed using (5.5e). 

The potential Q 2 d(z) is shown in Figure 5.1 Note that Q 2 d{z) is horizontal at 
z = 0, and monotonically decreasing in \z\. The negative of the slope \ = ~Q'id 
reaches a maximum of 



Xr, 



0.93305 at 0.59505. 



(5.8) 



The quantity Xmax plays a key role in our analysis of minimizers below. 

The Fourier transform of Q 2 d{z) can be evaluated exactly using the integral 
definition (5.5c) and interchanging the order of integration of z and w to obtain 

2tt 



Q2 D {k) = 



(1 + fc 2 ) 



2^/2 ' 



(5.9) 



which we note is positive, so local minimizers are global minimizers per the discussion 
in Section l2~3l 
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Fig. 5.1. Quasi-two-dimensional Laplace potential Q2D given by Q2D * s horizontal at 

■ 0, and monotonically decreasing in \z\. There are inflection points at ±z* given by ![5.8ty. 



5.2. Gravitational potential on the half-line with the quasi-2D-Laplace 
potential. We model a quasi-two-dimensional biological swarm with repulsive so- 
cial interactions of Laplace type and subject to an exogenous gravitational potential, 
F(x) — gx. The spatial coordinate x > describes the elevation above ground. 
Consequently, VL is the semi- infinite interval < x < 00. 



From Section 3.4 recall that for the one-dimensional model, 



p(x) = MS(x), 



(5.10) 



is a minimizer for some M, corresponding to all swarm members pinned by gravity to 
the ground. We consider this same solution as a candidate minimizer for the quasi- 
two-dimensional problem. In this case, p above is actually a minimizer for any mass 
M. To see this, we can compute A(x), 



A(x) 



Q(x - y)p{y) dy + F(x) = MQ 2 d(x) + gx. 



(5.11) 



Since A'(0) = g > 0, A increases away from the origin and hence p is at least a swarm 
minimizer. 

In fact, if M < Mi = g/xmaxi p(%) is a global minimizer because 



A'{x) = MQ' 2D (x) +g> g- M Xmax > 0, 



(5.12) 



which guarantees that A(x) is strictly increasing for x > as shown in Figure [5~2) [a) . 
Because it is strictly increasing, A(a;) > A(0) = A for x > 0. Given this fact, and 
additionally, since W2 > as previously shown, p is a global minimizer. This means 
that if an infinitesimal amount of mass is added anywhere in the system, it will descend 
to the origin. Consequently, we believe this solution is the global attractor (though 
we have not proven this). 

Note that while the condition M < Mi is sufficient for p to be a global mini- 
mizer, it is not necessary. As alluded above, it is not necessary that A(x) be strictly 
increasing, only that A(x) > A(0) for x > 0. This is the case for for M < Mi = g/x2, 
where \i ~ 0.79870. Figure [Ki^b) shows a case when Mi < M < M2. Although 
A(x) > A(0) for x > 0, A(x) has a local minimum. In this situation, although the so- 
lution with the mass concentrated at the origin is a global minimizer, it is not a global 
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attractor. We will see that a small amount of mass added near the local minimum of 
A(x) will create a swarm minimizer, which is dynamically stable to perturbations. 

Figure [5~2| c) shows the critical case when M — M%. In this case the local min- 
imum of A(x) at x — X2 ~ 1.11436 satisfies A(x2) = A(0) and A'(x2) = 0. Figure 
|5.2[ d) shows the case when M > M 2 and now A(x) < A(0) in the neighborhood of 
the minimum. In this case the solution with the mass concentrated at the origin is 
only a swarm minimizer; the energy of the system can be reduced by transporting 
some of the mass at the origin to the neighborhood of the local minimum. 

When M > M\ it is possible to construct a continuum of swarm minimizers. 
We have conducted a range of simulations for varying M and have measured two 
basic properties of the solutions. We set g = 1 and use N — 200 in all simulations 
of the discrete system. Initially, all the swarm members are high above the ground 
and we evolve the simulation to equilibrium. Figure [5T3] ^ a) measures the mass on the 
ground as a percentage of the total swarm mass. The horizontal blue line indicates 
(schematically) that for M < Mi, the equilibrium consists of all mass concentrated 
at the origin; as discussed above, this state is the global minimizer and (we believe) 
the global attractor. As mass is increased through Mi, the equilibrium is a swarm 
minimizer consisting of a classical swarm in the air separated from the origin, and 
some mass concentrated on the ground. As M increases, the proportion of mass 
located on the ground decreases monotonically. Figure |5.3[ b) visualizes the support 
of the airborne swarm, which exists only for M > Mi; the lower and upper data 
represent the coordinates of the bottom and top of the swarm, respectively. As mass 
is increased, the span of the swarm increases monotonically. 

As established above, when M > M\, swarm minimizers exist with two compo- 
nents. In fact, there is a continuum of swarm minimizers with different proportions 
of mass in the air and on the ground. Which minimizer is obtained in simulation 
depends on initial conditions. Figure 5.4 shows two such minimizers for g — 1 and 
M = 15 > M%, and the associated values of A(a;) (each obtained from a different 
initial condition). Recalling that for a swarm minimizer, each connected component 
of the swarm, A is constant, we define A(x) = A(0) = Ao for the grounded component 
and A (a;) = Ai for the airborne component. In Figure [BT^ ab) , 29.5% of the mass 
is contained in the grounded component. In this case, Ao < Ai indicating that the 
total energy could be reduced by transporting swarm members from the air to the 
ground. In contrast, in Figure [5T4| cd), 32.5% of the mass is contained in the grounded 
component. In this case, Ai < Ao indicating that the total energy could be reduced by 
transporting swarm members from the ground to the air. Note that by continuity, we 
believe a state exists where Ao = Ax, which would correspond to a global minimizer. 
However, this state is clearly not a global attractor and hence will not necessarily be 
achieved in simulation. 

We've demonstrated that for M > Mi one can construct a continuum of swarm 
minimizers with a gap between grounded and airborne components, and that for 
M > M-2 these solutions can have a lower energy than the state with the density 
concentrated solely on the ground. By contrast with the one-dimensional system of 
Section |3.4| in which no gap is observed, these gap states appear to be the generic 
configuration for sufficiently large mass in the quasi-two-dimensional system. We 
conclude that dimensionality is crucial element for the formation of the bubble-like 
shape of real locust swarms. 



6. Conclusions. In this paper we deeveloped a framework for studying equilib- 
rium solutions for swarming problems. We related the discrete swarming problem to 
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Fig. 5.2. A(a:) in \2. 24\ ) for the half-line f! = [0, oo] with exogenous gravitational potential 
F(x) = gx with g = 1, and endogenous interactions given by the quasi-two-dimensional Laplace 
potential {5..^ . (a) M = 1 < Mi, in which case the S- concentration at the x = (represented by 
the dot) is a global miniraizer and (we believe) a global attractor. (b) M2 > M = 1.15 > Mi, in 
which case the 5 -concentration is a global miniraizer but not a global attractor. Another equilibrium 
solution can be found by transferring a small amount of mass from the origin to the local minimum 
in A(a;). (c) M = M2 ~ 1.25204. This is the boundary case; for M > M2, the 5 -concentration at the 
origin is no longer a global minimizer, although it remains a swarm minimizer. (d) M = 1.4 > M2, 
in which case the 5 -concentration at the origin is now only a swarm minimizer. 



an associated continuum model. This continuum model has an energy formulation 
which enables analysis equilibrium solutions and their stability. We derived condi- 
tions for an equilibrium solution to be a local minimizer, a global minimizer, and/or 
a swarm minimizer, that is, stable to infinitesimal Lagrangian deformations of the 
mass. 

We found many examples of compactly supported equilibrium solutions, which 
may be discontinuous at the boundary of the support. In addition, when a bound- 
ary of the support coincides with the domain boundary, a minimizer may contain a 
^-concentration there. For the case of exogenous repulsion modeled by the Laplace 
potential, we computed three example equilibria. On a bounded domain, the mini- 
mizer is a constant density profile with (^-functions at each end. On a half-line with 
an exogenous gravitational potential, the minimizer is a compactly supported lin- 
ear density profile with a <5-function at the origin. In free space with an exogenous 
quadratic potential, the minimizer is a compactly supported inverted parabola with 
jump discontinuities at the endpoints. Each of the aforementioned solutions is also a 
global minimizer. 
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Fig. 5.3. Numerical simulations of \2. Ity for the half-line Q = [0, oo] with exogenous gravita- 
tional pot entia l F(x) = gx and endogenous interactions given by the quasi-two-dimensional Laplace 
potential \5.4\ - (a) Mass on the ground as a percentage of the total swarm mass. The horizontal blue 
line indicates (schematically) that for M < Mi, the equilibrium consists of all mass concentrated at 
the origin. As mass is increased through Mi, a classical swarm exists in the air separated from the 
origin, and the proportion of mass located on the ground decreases monotonically. (b ) Support of 
the airborne swarm, which exists only for M > Mi; the lower and upper data represent the coordi- 
nates of the bottom and top of the swarm, respectively. The bottom of the swarm is located a finite 
distance from the origin, which is represented as the horizontal dotted line. As mass is increased, 
the span of the swarm increases monotonically, with the location of the bottom remaining constant 
within the error of our Lagrangian numerical approximation. For all computations, we take g = 1 
and N = 200 swarm members in the discrete simulation. We use a linear spacing in M for larger 
values of M , but a logarithmic spacing for values close to Mi in order to resolve the bifurcation. 



To extend the results above, we also found analytical solutions for exogenous 
attractive-repulsive forces, modeled with the Morse potential. In the case that the so- 
cial force was in the catastrophic statistical mechanical regime, we found a compactly 
supported solution whose support is independent of the total population mass. This 
means that within the modeling assumptions, swarms become denser with increasing 
mass. For the case of an H-stable social force, there is no equilibrium solution on an 
infinite domain. On a finite domain, mass is partitioned between a classical solution 
in the interior and ^-concentrations on the boundary. 

We recall that for the locust model of [33] (see Figure |1.1[ ) a concentration of lo- 
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Fig. 5.4. Two example swarm minimizers for the quasi-two-dimensional Laplace potential 
j"5.4P on the half-line fl = [0, oo] with exogenous gravitational potential F(x) = gx. We take g = 1 
and N = 200 swarm members for our simulations, but choose different initial conditions. Each 
minimizer consists of two swarm components, namely a concentration of mass at the origin and a 
finite-sized swarm separated from the origin by a gap. (ab) 29.5% of the mass is in the grounded 
component, (a) Shows A(x). The dot and cross-hatched bar schematically indicate the support of the 
two components. Since Ao = A(0) < Ai where Ai is A(x) in the airborne component, we conclude 
that the total energy could be reduced by transporting mass from the air to the ground, (b) The 
swarm density p(x). The "lollipop" represents a superposition of 59 swarm members, having total 
mass 59 • 1/200. Density within the airborne component is given only at the Lagrangian gridpoints. 
The horizontal arrow helps demarcate its support, (cd) Like (ab), but here 32.5% of the mass is 
contained in the grounded group. Since Ao < Ai, the total energy could be reduced by transporting 
mass from the ground to the air. Note that although both states could have their energy reduced 
by transporting mass between components, each one is still a dynamically stable equilibrium. We 
believe that these are just two representative examples of a continuum of dynamically stable swarm 
equilibria all with the same total mass. 



custs occurs on the ground, with a seemingly classical component above, separated by 
a gap. None of the one-dimensional solutions (for the Laplace and Morse potentials) 
discussed above contain a gap, that is, multiple swarm components that are spatially 
disconnected, suggesting that this configuration is intrinsically two-dimensional. To 
study this configuration, we computed a quasi-two-dimensional potential correspond- 
ing to a horizontally uniform swarm. We demonstrated numerically that for a wide 
range of parameters, there exists a continuous family of swarm minimizers that consist 
of a concentration on the ground and a disconnected, classical component in the air, 
reminiscent of our earlier numerical studies of a discrete locust swarm model. 

We believe that the analytical solutions we found provide a sampling of the rich 
tapestry of equilibrium solutions that manifest in the general model we have consid- 
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ered, and in nature. We hope that these solutions will inspire further analysis and 
guide future modeling efforts. 
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